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Abstract 

Sparse data models, where data is assumed to be well represented as a linear combination of a few elements 
from a dictionary, have gained considerable attention in recent years, and their use has led to state-of-the-art results 
in many signal and image processing tasks. It is now well understood that the choice of the sparsity regularization 
term is critical in the success of such models. Based on a codelength minimization interpretation of sparse coding, 
and using tools from universal coding theory, we propose a framework for designing sparsity regularization terms 
which have theoretical and practical advantages when compared to the more standard £o or i\ ones. The presentation 
of the framework and theoretical foundations is complemented with examples that show its practical advantages in 
image denoising, zooming and classification. 

I. Introduction 

Sparse modeling calls for constructing a succinct representation of some data as a combination of a few typical 
patterns {atoms) learned from the data itself. Significant contributions to the theory and practice of learning such 
collections of atoms (usually called dictionaries or codebooks), e.g., [1], [14], [33], and of representing the actual 
data in terms of them, e.g., [8], [11], [12], have been developed in recent years, leading to state-of-the-art results in 
many signal and image processing tasks [24], [26], [27], [34]. We refer the reader for example to [ ] for a recent 
review on the subject. 

A critical component of sparse modeling is the actual sparsity of the representation, which is controlled by a 
regularization term (regularizer for short) and its associated parameters. The choice of the functional form of the 
regularizer and its parameters is a challenging task. Several solutions to this problem have been proposed in the 
literature, ranging from the automatic tuning of the parameters [20] to Bayesian models, where these parameters 
are themselves considered as random variables [17], [20], [51]. In this work we adopt the interpretation of sparse 
coding as a codelength minimization problem. This is a natural and objective method for assessing the quality 
of a statistical model for describing given data, and which is based on the Minimum Description Length (MDL) 
principle [37]. In this framework, the regularization term in the sparse coding formulation is interpreted as the cost 
in bits of describing the sparse linear coefficients used to reconstruct the data. Several works on image coding 
using this approach were developed in the 1990's under the name of "complexity -based" or "compression-based" 
coding, following the popularization of MDL as a powerful statistical modeling tool [9], [31], [40]. The focus on 
these early works was in denoising using wavelet basis, using either generic asymptotic results from MDL or fixed 
probability models, in order to compute the description length of the coefficients. A later, major breakthrough in 
MDL theory was the adoption of universal coding tools to compute optimal codelengths. In this work, we improve 
and extend on previous results in this line of work by designing regularization terms based on such universal codes 
for image coefficients, meaning that the codelengths obtained when encoding the coefficients of any (natural) image 
with such codes will be close to the shortest codelengths that can be obtained with any model fitted specifically for 
that particular instance of coefficients. The resulting framework not only formalizes sparse coding from the MDL 
and universal coding perspectives but also leads to a family of universal regularizes which we show to consistently 
improve results in image processing tasks such as denoising and classification. These models also enjoy several 
desirable theoretical and practical properties such as statistical consistency (in certain cases), improved robustness 
to outliers in the data, and improved sparse signal recovery (e.g., decoding of sparse signals from a compressive 
sensing point of view [5]) when compared with the traditional £o and i\ -based techniques in practice. These models 
also yield to the use of a simple and efficient optimization technique for solving the corresponding sparse coding 
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problems as a series of weighted l\ subproblems, which in turn, can be solved with off-the-shelf algorithms such 
as LARS [ ] or 1ST [11]. Details are given in the sequel. 

Finally, we apply our universal regularizers not only for coding using fixed dictionaries, but also for learning the 
dictionaries themselves, leading to further improvements in all the aforementioned tasks. 

The remainder of this paper is organized as follows: in Section II we introduce the standard framework of sparse 
modeling. Section III is dedicated to the derivation of our proposed universal sparse modeling framework, while 
Section IV deals with its implementation. Section V presents experimental results showing the practical benefits 
of the proposed framework in image denoising, zooming and classification tasks. Concluding remarks are given in 
Section VI. 

II. Sparse modeling and the need for better models 
Let X <E R MxN be a set of N column data samples Xj <E R M , D E ^MxK a dictionary of K column atoms 
dk E M M , and A E R Kx7V , aj E a set of reconstruction coefficients such that X = D A. We use a^ to denote 
the k-th row of A, the coefficients associated to the k-th atom in D. For each j = 1, . . . , N we define the active 
set of a.j as Aj = {k : a^j ^ 0, 1 < k < K}, and ||aj|| = \Aj\ as its cardinality. The goal of sparse modeling 
is to design a dictionary D such that for all or most data samples Xj, there exists a coefficients vector aj such 
that Xj « Daj and ||aj|| is small (usually below some threshold L <C K). Formally, we would like to solve the 
following problem 

N 

min ^^^(aj) s.t. ||xj — D aj ||^ < e, j = 1, . . . , JV, (1) 

where ip(-) is a regularization term which induces sparsity in the columns of the solution A. Usually the constraint 
||dfc|| 2 < 1, k = is added, since otherwise we can always decrease the cost function arbitrarily by 

multiplying D by a large constant and dividing A by the same constant. When D is fixed, the problem of finding 
a sparse aj for each sample Xj is called sparse coding, 

a.j = argmin V>( a j) s -t- ll x j — Da^H^ < e. (2) 

Among possible choices of ^( ) are the £o pseudo-norm, ^( ) = || || , and the i\ norm. The former tries to 
solve directly for the sparsest aj, but since it is non-convex, it is commonly replaced by the i\ norm, which is 
its closest convex approximation. Furthermore, under certain conditions on (fixed) D and the sparsity of a j 9 the 
solutions to the £o and ^i-based sparse coding problems coincide (see for example [ ]). The problem (1) is also 
usually formulated in Lagrangian form, 

N 

along with its respective sparse coding problem when D is fixed, 

aj = argmin ||xj -Da^i A^(a). (4) 

Even when the regularizer is convex, the sparse modeling problem, in any of its forms, is jointly non-convex in 
(D, A). Therefore, the standard approach to find an approximate solution is to use alternate minimization: starting 
with an initial dictionary T)(°\ we minimize (3) alternatively in A via (2) or (4) (sparse coding step), and D 
(dictionary update step). The sparse coding step can be solved efficiently when ?/>(•) = H*^ using for example 1ST 
[11] or LARS [12], or with OMP [28] when = || || . The dictionary update step can be done using for example 
MOD [14] or K-SVD [ ]. 

A. Interpretations of the sparse coding problem 

We now turn our attention to the sparse coding problem: given a fixed dictionary D, for each sample vector Xj, 
compute the sparsest vector of coefficients aj that yields a good approximation of Xj. The sparse coding problem 
admits several interpretations. What follows is a summary of these interpretations and the insights that they provide 
into the properties of the sparse models that are relevant to our derivation. 
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1) Model selection in statistics: Using the £q norm as xjj(') in (4) is known in the statistics community as the 
Akaike's Information Criterion (AIC) when A = 1, or the Bayes Information Criterion (BIC) when A = ^logM, 
two popular forms of model selection (see [ , Chapter 7]). In this context, the t\ regularizer was introduced in 
[43], again as a convex approximation of the above model selection methods, and is commonly known (either in 
its constrained or Lagrangian forms) as the Lasso. Note however that, in the regression interpretation of (4), the 
role of D and X is very different. 

2) Maximum a posteriori: Another interpretation of (4) is that of a maximum a posteriori (MAP) estimation of 
Rj in the logarithmic scale, that is 

a.j = argmax {log P(a|xj)} = argmax {logP(xj|a) + logP(a)} 

= argmin {— logP(xj|a) — logP(a)}, (5) 

where the observed samples Xj are assumed to be contaminated with additive, zero mean, IID Gaussian noise with 
variance a 2 , P(xj|a) oc " Xj_Da " 2 , and a prior probability model on a with the form P(a) oc e~°^^ is 

considered. The energy term in Equation (4) follows by plugging the previous two probability models into (5) and 
factorizing 2a 2 into A = 2a 2 9. According to (5), the i\ regularizer corresponds to an IID Laplacian prior with 
mean and inverse-scale parameter 0, P(a) = fl^Li 0e~ e ^ = 9 K e~°^ a ^ 9 which has a special meaning in signal 
processing tasks such as image or audio compression. This is due to the widely accepted fact that representation 
coefficients derived from predictive coding of continuous-valued signals, and, more generally, responses from zero- 
mean filters, are well modeled using Laplacian distributions. For example, for the special case of DCT coefficients 
of image patches, an analytical study of this phenomenon is provided in [25], along with further references on the 
subject. 

3) Codelength minimization: Sparse coding, in all its forms, has yet another important interpretation. Suppose 
that we have a fixed dictionary D and that we want to use it to compress an image, either losslessly by encoding the 
reconstruction coefficients A and the residual X — DA, or in a lossy manner, by obtaining a good approximation 
X & DA and encoding only A. Consider for example the latter case. Most modern compression schemes consist 
of two parts: a probability assignment stage, where the data, in this case A, is assigned a probability P(A), and 
an encoding stage, where a code C(A) of length L(A) bits is assigned to the data given its probability, so that 
L(A) is as short as possible. The techniques known as Arithmetic and Huffman coding provide the best possible 
solution for the encoding step, which is to approximate the Shannon ideal codelength L(A) = — logP(A) [10, 
Chapter 5]. Therefore, modern compression theory deals with finding the coefficients A that maximize P(A), or, 
equivalently, that minimize — logP(A). Now, to encode X lossily, we obtain coefficients A such that each data 
sample Xj is approximated up to a certain £2 distortion e, ||xj — Da^H^ < e. Therefore, given a model P(a) for a 
vector of reconstruction coefficients, and assuming that we encode each sample independently, the optimum vector 
of coefficients rj for each sample Xj will be the solution to the optimization problem 

aj = argmin — logP(a) s.t. ||xj — Da^H^ < e, (6) 

which, for the choice P(a) oc e~^^ coincides with the error constrained sparse coding problem (2). Suppose now 
that we want lossless compression. In this case we also need to encode the reconstruction residual Xj — Daj. Since 
P(x, a) = P(x|a)P(a), the combined codelength will be 

L(x j ,a j ) = -logP(x j ,a j ) = - logP(x j |a j ) -logP(a j ). (7) 

Therefore, obtaining the best coefficients rj amounts to solving min a L(xj, a^), which is precisely the MAP 
formulation of (5), which in turn, for proper choices of P(x|a) and P(a), leads to the Lagrangian form of sparse 
coding (4). 1 

laplacian models, as well as Gaussian models, are probability distributions over R, characterized by continuous probability density 
functions, f(a) = F'(a), F(a) = P(x < a). If the reconstruction coefficients are considered real numbers, under any of these distributions, 
any instance of A G R KxiV will have measure 0, that is, P(A) = 0. In order to use such distributions as our models for the data, 
we assume that the coefficients in A are quantized to a precision A, small enough for the density function f(a) to be approximately 
constant in any interval [a — A/2, a + A/2], a G R, so that we can approximate P(a) ~ A/ (a), a G R. Under these assumptions, 
— logP(a) ~ — log f(a) — log A, and the effect of A on the codelength produced by any model is the same. Therefore, we will omit A 
in the sequel, and treat density functions and probability distributions interchangeably as P(-). Of course, in real compression applications, 
A needs to be tuned. 
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Fig. 1. Standard 8x8 DCT dictionary (a), global empirical distribution of the coefficients in A (b, log scale), empirical distributions of the 
coefficients associated to each of the K = 64 DCT atoms (c, log scale). The distributions in (c) have a similar heavy tailed shape (heavier 
than Laplacian), but the variance in each case can be significantly different, (d) Histogram of the K = 64 different 9k values obtained by 
fitting a Laplacian distribution to each row of A. Note that there are significant occurrences between 9 — 5 to 9 — 25. The coefficients A 
used in (b-d) were obtained from encoding 10 6 8x8 patches (after removing their DC component) randomly sampled from the Pascal 2006 
dataset of natural images [ ]. (e) Histograms showing the spatial variability of the best local estimations of 9k for a few rows of A across 
different regions of an image. In this case, the coefficients A correspond to the sparse encoding of all 8x8 patches from a single image, in 
scan-line order. For each k, each value of 9k was computed from a random contiguous block of 250 samples from . The procedure was 
repeated 4000 times to obtain an empirical distribution. The wide supports of the empirical distributions indicate that the estimated 9 can 
have very different values, even for the same atom, depending on the region of the data from where the coefficients are taken. 



As one can see, the codelength interpretation of sparse coding is able to unify and interpret both the constrained 
and unconstrained formulations into one consistent framework. Furthermore, this framework offers a natural and 
objective measure for comparing the quality of different models P(x|a) and P(a) in terms of the codelengths 
obtained. 

4) Remarks on related work: As mentioned in the introduction, the codelength interpretation of signal coding 
was already studied in the context of orthogonal wavelet-based denoising. An early example of this line of work 
considers a regularization term which uses the Shannon Entropy function J2Pi^°&Pi t0 gi ye a measure of the 
sparsity of the solution [9]. However, the Entropy function is not used as measure of the ideal codelength for 
describing the coefficients, but as a measure of the sparsity (actually, group sparsity) of the solution. The MDL 
principle was applied to the signal estimation problem in [ ]. In this case, the codelength term includes the 
description of both the location and the magnitude of the nonzero coefficients. Although a pioneering effort, the 
model assumed in [40] for the coefficient magnitude is a uniform distribution on [0, 1], which does not exploit a 
priori knowledge of image coefficient statistics, and the description of the support is slightly wasteful. Furthermore, 
the codelength expression used is an asymptotic result, actually equivalent to BIC (see Section II-A1) which can 
be misleading when working with small sample sizes, such as when encoding small image patches, as in current 
state of the art image processing applications. The uniform distribution was later replaced by the universal code 
for integers [ ] in [ ]. However, as in [ ], the model is so general that it does not perform well for the specific 
case of coefficients arising from image decompositions, leading to poor results. In contrast, our models are derived 
following a careful analysis of image coefficient statistics. Finally, probability models suitable to image coefficient 
statistics of the form P(a) oc e _ ' a ^ (known as generalized Gaussians) were applied to the MDL-based signal 
coding and estimation framework in [ ]. The justification for such models is based on the empirical observation 
that sparse coefficients statistics exhibit "heavy tails" (see next section). However, the choice is ad hoc and no 
optimality criterion is available to compare it with other possibilities. Moreover, there is no closed form solution 
for performing parameter estimation on such family of models, requiring numerical optimization techniques. In 
Section III, we derive a number of probability models for which parameter estimation can be computed efficiently 
in closed form, and which are guaranteed to optimally describe image coefficients. 

B. The need for a better model 

As explained in the previous subsection, the use of the t\ regularizer implies that all the coefficients in A 
share the same Laplacian parameter 9. However, as noted in [ ] and references therein, the empirical variance of 
coefficients associated to different atoms, that is, of the different rows of A, varies greatly with k = 1 . . . , K. 
This is clearly seen in Figures l(a-c), which show the empirical distribution of DCT coefficients of 8x8 patches. 
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As the variance of a Laplacian is 2/9 2 , different variances indicate different underlying 9. The histogram of the set 
j^, k = 1, . . . , if j of estimated Laplacian parameters for each row k, Figure 1(d), shows that this is indeed the 

case, with significant occurrences of values of 9 in a range of 5 to 25. 

The straightforward modification suggested by this phenomenon is to use a model where each row of A has 
its own weight associated to it, leading to a weighted i\ regularizes However, from a modeling perspective, this 
results in K parameters to be adjusted instead of just one, which often results in poor generalization properties. For 
example, in the cases studied in Section V, even with thousands of images for learning these parameters, the results 
of applying the learned model to new images were always significantly worse (over ldB in estimation problems) 
when compared to those obtained using simpler models such as an unweighted i\. 2 One reason for this failure may 
be that real images, as well as other types of signals such as audio samples, are far from stationary. In this case, 
even if each atom k is associated to its own 9^ (A&), the optimal value of 9^ can have significant local variations 
at different positions or times. This effect is shown in Figure 1(e), where, for each k, each 9^ was re-estimated 
several times using samples from different regions of an image, and the histogram of the different estimated values 
of Ok was computed. Here again we used the DCT basis as the dictionary D. 

The need for a flexible model which at the same time has a small number of parameters leads naturally to Bayesian 
formulations where the different possible are "marginalized out" by imposing an hyper-prior distribution on A, 
sampling A using its posterior distribution, and then averaging the estimates obtained with the sampled sparse- 
coding problems. Examples of this recent line of work, and the closely related Bayesian Compressive Sensing, are 
developed for example in [23], [44], [49], [48]. Despite of its promising results, the Bayesian approach is often 
criticized due to the potentially expensive sampling process (something which can be reduced for certain choices 
of the priors involved [23]), arbitrariness in the choice of the priors, and lack of proper theoretical justification for 
the proposed models [48]. 

In this work we pursue the same goal of deriving a more flexible and accurate sparse model than the traditional 
ones, while avoiding an increase in the number of parameters and the burden of possibly solving several sampled 
instances of the sparse coding problem. For this, we deploy tools from the very successful information-theoretic 
field of universal coding, which is an extension of the compression scenario summarized above in Section II-A, 
when the probability model for the data to be described is itself unknown and has to be described as well. 

III. Universal models for sparse coding 

Following the discussion in the preceding section, we now have several possible scenarios to deal with. First, 
we may still want to consider a single value of 9 to work well for all the coefficients in A, and try to design a 
sparse coding scheme that does not depend on prior knowledge on the value of 9. Secondly, we can consider an 
independent (but not identically distributed) Laplacian model where the underlying parameter 9 can be different 
for each atom k = 1, . . . , K. In the most extreme scenario, we can consider each single coefficient a^j in A 
to have its own unknown underlying 0jy and yet, we would like to encode each of these coefficients (almost) as if 
we knew its hidden parameter. 

The first two scenarios are the ones which fit the original purpose of universal coding theory [ ], which is 
the design of optimal codes for data whose probability models are unknown, and the models themselves are to be 
encoded as well in the compressed representation. 

Now we develop the basic ideas and techniques of universal coding applied to the first scenario, where the 
problem is to describe A as an IID Laplacian with unknown parameter 9. Assuming a known parametric form for 
the prior, with unknown parameter 9, leads to the concept of a model class. In our case, we consider the class 
M = {P(A\9) : 9 G 0} of all iid Laplacian models over A e R KxN , where 

N K 

p{M°) = II II p ^kj\e) = 9e-°^ 

and 6 C R + . The goal of universal coding is to find a probability model Q(A) which can fit A as well as the 
model in M that best fits A after having observed it. A model Q(A) with this property is called universal (with 
respect to the model M). 

2 Note that this is the case when the weights are found by maximum likelihood. Other applications of weighted £i regularizers, using other 
types of weighting strategies, are known to improve over £± -based ones for certain applications (see e.g. [51]). 
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For simplicity, in the following discussion we consider the coefficient matrix A to be arranged as a single long 
column vector of length n = K xN, a = (ai, . . . , a n ). We also use the letter a without sub-index to denote the 
value of a random variable representing coefficient values. 

First we need to define a criterion for comparing the fitting quality of different models. In universal coding theory 
this is done in terms of the codelengths L(a) required by each model to describe a. 

If the model consists of a single probability distribution P(-), we know from Section II- A3 that the optimum 
codelength corresponds to Lp(a) = — logP(a). Moreover, this relationship defines a one-to-one correspondence 
between distributions and codelengths, so that for any coding scheme Lg(a), Q(a) = 2~ Lq ^\ Now suppose 
that we are restricted to a class of models M, and that we need choose the model P E M that assigns the 
shortest codelength to a particular instance of a. We then have that P is the model in M that assigns the maximum 
probability to a. For a class M parametrized by 9, this corresponds to P = P(a|#(a)), where 0(a) is the maximum 
likelihood estimator (MLE) of the model class parameter 9 given a (we will usually omit the argument and just 
write 9). Unfortunately, we also need to include the value of 9 in the description of a for the decoder to be able to 
reconstruct it from the code C(a). Thus, we have that any model Q(a) inducing valid codelengths Lq(sl) will have 
Lq(sl) > — logP(a|0). The overhead of Lq(s.) with respect to — logP(a|0) is known as the codelength regret, 

U( a ,Q) := L Q (a) - (-logP(a|0(a))) = -logQ(a) + logP(a|0(a))). 

A model Q( a ) ( or > more precisely, a sequence of models, one for each data length n) is called universal if 7£(a, Q) 
grows sublinearly in n for all possible realizations of a, that is ^7£(a, Q) — > , Va E M n , so that the codelength 
regret with respect to the MLE becomes asymptotically negligible. 

There are a number of ways to construct universal probability models. The simplest one is the so called two- 
part code, where the data is described in two parts. The first part describes the optimal parameter 0(a) and the 
second part describes the data according to the model with the value of the estimated parameter 9, P(a|0(a)). For 
uncountable parameter spaces 6, such as a compact subset of R, the value of 9 has to be quantized in order to be 
described with a finite number of bits d. We call the quantized parameter 9^ The regret for this model is thus 

ft(a, Q) = L{9 d ) + L(a|0 d ) - L(a|0) = L{6 d ) - logP(a|0 d ) - (-logP(a|0)). 

The key for this model to be universal is in the choice of the quantization step for the parameter 9, so that both 
its description L(0d), and the difference — logP(a|0^) — (— logP(a|0)), grow sublinearly. This can be achieved 
by letting the quantization step shrink as 0(l/y/ri) [ ], thus requiring d = 0(0.5 log n) bits to describe each 
dimension of 9^ This gives us a total regret for two-part codes which grows as d ™( e ) logn, where dim(0) is the 
dimension of the parameter space 6. 

Another important universal code is the so called Normalized Maximum Likelihood (NML) [ ]. In this case the 
universal model Q*(a) corresponds to the model that minimizes the worst case regret, 

Q*(a) = minmax{— logQ(a) + log P(a|0(a))}, 
Q a 

which can be written in closed form as Q*(a) = P q^^ > where the normalization constant 

C(M,n) := J2 ^(a|0(a))da 

is the value of the minimax regret and depends only on M and the length of the data n? Note that the NML model 
requires C(M, n) to be finite, something which is often not the case. 

The two previous examples are good for assigning a probability to coefficients a that have already been computed, 
but they cannot be used as a model for computing the coefficients themselves since they depend on having observed 
them in the first place. For this and other reasons that will become clearer later, we concentrate our work on a third 
important family of universal codes derived from the so called mixture models (also called Bayesian mixtures). In 

3 The minimax optimality of Q* (a) derives from the fact that it defines a complete uniquely decodable code for all data a of length n, that 
is, it satisfies the Kraft inequality with equality. ^ aeRri 2~ L Q*^ — 1. Since every uniquely decodable code with lengths {Lq{sl) : a £ R n } 
must satisfy the Kraft inequality (see [ , Chapter 5]), if there exists a value of a such that Lq(a) < Lq* (a) (that is 2~ L Q (a) > 2~ L Q* (a) ), 
then there exists a vector a 7 for which Lq(8l) > Lq* (a ; ) for the Kraft inequality to hold. Therefore the regret of Q for a ; is necessarily 
greater than C(Ai,n), which shows that Q* is minimax optimal. 
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a mixture model, Q(a) is a convex mixture of all the models P(sl\0) in M, indexed by the model parameter 9, 
Q(a) = J e P(a\0)w(0)dO, where w{6) specifies the weight of each model. Being a convex mixture implies that 
w{9) > and J e w(9)d9 = 1, thus w{9) is itself a probability measure over 6. We will restrict ourselves to the 
particular case when a is considered a sequence of independent random variables, 4 

Q(a) = I]Q j (a j ) > Qjifij) = / ^(a#H(0)^, (8) 

where the mixing function wj(6) can be different for each sample j. An important particular case of this scheme 
is the so called Sequential Bayes code, in which wj(9) is computed sequentially as a posterior distribution based 
on previously observed samples, that is Wj{9) = P{9\a\ 1 . . . , a n _i) [ , Chapter 6]. In this work, for simplicity, 
we restrict ourselves to the case where Wj(6) = w(6) is the same for all j. The result is an IID model where the 
probability of each sample aj is a mixture of some probability measure over R, 

Qj{a,j) = Q{a,j) = [ P( aj \0)w(0)(W,\/j = l,...,N. (9) 

A well known result for IID mixture (Bayesian) codes states that their asymptotic regret is 0( dm ^ e ) logn), thus 
stating their universality, as long as the weighting function w{6) is positive, continuous and unimodal over 6 (see 
for example [21, Theorem 8.1],[ ]). This gives us great flexibility on the choice of a weighting function w{6) 
that guarantees universality. Of course, the results are asymptotic and the o(log n) terms can be large, so that the 
choice of w{9) can have practical impact for small sample sizes. 

In the following discussion we derive several IID mixture models for the Laplacian model class M. For this 
purpose, it will be convenient to consider the corresponding one-sided counterpart of the Laplacian, which is the 
exponential distribution over the absolute value of the coefficients, \a\, and then symmetrize back to obtain the final 
distribution over the signed coefficients a. 



A. The conjugate prior 

In general, (9) can be computed in closed form if w(6) is the conjugate prior of P(a\9). When P{a\9) is an 
exponential (one-sided Laplacian), the conjugate prior is the Gamma distribution, 

w (9\k, /3) = T{K)- 1 9 K - l ^e-^ e , 9 e R + , 

where n and (3 are its shape and scale parameters respectively. Plugging this in (9) we obtain the Mixture of 
exponentials model (moe), which has the following form (see Appendix A for the full derivation), 

Q M0E (a|/3, k) = K(3*(a + a e M + . (10) 

With some abuse of notation, we will also denote the symmetric distribution on a as MOE, 

Q M0E (a|/3, «) = \^{\a\ + /3)"^ +1 ), a e R. (11) 

Although the resulting prior has two parameters to deal with instead of one, we know from universal coding 
theory that, in principle, any choice of k and f3 will give us a model whose codelength regret is asymptotically 
small. 

Furthermore, being IID models, each coefficient of a itself is modeled as a mixture of exponentials, which makes 
the resulting model over a very well suited to the most flexible scenario where the "underlying" 9 can be different 
for each aj. In Section V-B we will show that a single MOE distribution can fit each of the K rows of A better than 
K separate Laplacian distributions fine-tuned to these rows, with a total of K parameters to be estimated. Thus, 
not only we can deal with one single unknown 9, but we can actually achieve maximum flexibility with only two 
parameters (k and (3). This property is particular of the mixture models, and does not apply to the other universal 
models presented. 

4 More sophisticated models which include dependencies between the elements of a are out of the scope of this work. 
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Finally, if desired, both n and (3 can be easily estimated using the method of moments (see Appendix A). Given 
sample estimates of the first and second non-central moments, jli = ^ Y^j=i \ a j\ an< ^ A 2 = \ Y^j=i l a j| 2 ' we have 
that 

k = 2 (#2 - Ai)/(A2 - 2/if) and /3 = (ft — l)/ii. (12) 

When the MOE prior is plugged into (5) instead of the standard Laplacian, the following new sparse coding 
formulation is obtained, 

K 

dj = argmin ||xj - Da^ + A M0E log (|a fe | + (3) , (13) 

fc=i 

where A M0E = 2a 2 (k, + 1). An example of the MOE regularizer, and the thresholding function it induces, is shown 
in Figure 2 (center column) for k = 2.5,(3 = 0.05. Smooth, differentiate non-convex regularizers such as the 
one in in (13) have become a mainstream robust alternative to the t\ norm in statistics [16], [51]. Furthermore, 
it has been shown that the use of such regularizers in regression leads to consistent estimators which are able to 
identify the relevant variables in a regression model (oracle property) [16]. This is not always the case for the i\ 
regularizer, as was proved in [ ]. The MOE regularizer has also been recently proposed in the context of compressive 
sensing [ ], where it is conjectured to be better than the ^i-term at recovering sparse signals in compressive sensing 
applications. 5 This conjecture was partially confirmed recently for non-convex regularizers of the form ^(a) = ||a|| r 
with < r < 1 in [39], [18], and for a more general family of non-convex regularizers including the one in (13) 
in [47]. In all cases, it was shown that the conditions on the sensing matrix (here D) can be significantly relaxed 
to guarantee exact recovery if non-convex regularizers are used instead of the i\ norm, provided that the exact 
solution to the non-convex optimization problem can be computed. In practice, this regularizer is being used with 
success in a number of applications here and in [7], [46]. 6 Our experimental results in Section V provide further 
evidence on the benefits of the use of non-convex regularizers, leading to a much improved recovery accuracy of 
sparse coefficients compared to i\ and £q. We also show in Section V that the MOE prior is much more accurate 
than the standard Laplacian to model the distribution of reconstruction coefficients drawn from a large database of 
image patches. We also show in Section V how these improvements lead to better results in applications such as 
image estimation and classification. 



B. The Jeffreys prior 

The Jeffreys prior for a parametric model class M = {P(a\0), 9 E O}, is defined as 



where is the determinant of the Fisher information matrix 



1(0) = E_ 



1 



|^log P(a\0) 



(15) 



The Jeffreys prior is well known in Bayesian theory due to three important properties: it virtually eliminates 
the hyper-parameters of the model, it is invariant to the original parametrization of the distribution, and it is a 
"non-informative prior," meaning that it represents well the lack of prior information on the unknown parameter 9 
[3]. It turns out that, for quite different reasons, the Jeffreys prior is also of paramount importance in the theory 
of universal coding. For instance, it has been shown in [2] that the worst case regret of the mixture code obtained 
using the Jeffreys prior approaches that of the NML as the number of samples n grows. Thus, by using Jeffreys, one 
can attain the minimum worst case regret asymptotically, while retaining the advantages of a mixture (not needing 
hindsight of a), which in our case means to be able to use it as a model for computing a via sparse coding. 

For the exponential distribution we have that 1(9) = p. Clearly, if we let 6 = (0, 00), the integral in (14) 
evaluates to 00. Therefore, in order to obtain a proper integral, we need to exclude and 00 from 6 (note that 

5 In [6], the logarithmic regularizer arises from approximating the £0 pseudo-norm as an t\ -normalized element- wise sum, without the 
insight and theoretical foundation here reported. 

6 While these works support the use of such non-convex regularizers, none of them formally derives them using the universal coding 
framework as in this paper. 
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a a a a a a 



Fig. 2. Left to right: £± (green), MOE (red) and JOE (blue) regularizers and their corresponding thresholding functions thres(x) := 
argmin a {(x — a) 2 + A^(|a|)}. The unbiasedness of MOE is due to the fact that large coefficients are not shrank by the thresholding 
function. Also, although the JOE regularizer is biased, the shrinkage of large coefficients can be much smaller than the one applied to small 
coefficients. 



this was not needed for the conjugate prior). We choose to define 6 = [0i,02]» < 9\ < 02 < oo, leading to 
The resulting mixture, after being symmetrized around 0, has the following form (see Appendix A): 

- ^m^)R (^'"" " *~ ,M ) • ° £ K+ < 16 > 

We refer to this prior as a Jeffreys mixture of exponentials (JOE), and again overload this acronym to refer to 
the symmetric case as well. Note that although Q JO e is not defined for a = 0, its limit when a — >> is finite and 
evaluates to 2 in(fl~ 2 /6> 1 ) - Thus, by defining Q JO e(0) = 2 in(^ 2 /6> 1 ) ? we °btain a prior that is well defined and continuous 
for all a £ R. When plugged into (5), we get the JOE-based sparse coding formulation, 

K 

mm || Xj - Da||* + A J0E ^{log \a k \ - log^ 1 ^ - e"^)}, (17) 

k=i 

where, according to the convention just defined for Q JO e(0), we define ^ JO e(0) := log(#2 — 9\). According to the 
MAP interpretation we have that A J0E = 2cr 2 , coming from the Gaussian assumption on the approximation error as 
explained in Section II-A. 

As with MOE, the JOE-based regularizer, t/'joeG) = — logQjo E (*)> is continuous and differentiable in R + , and its 
derivative converges to a finite value at zero, lim a ^o % E\ a ) = el-el ' ^ s we w ^ see ^ ater * n ^ ect i° n FV, these 
properties are important to guarantee the convergence of sparse coding algorithms using non-convex priors. Note 
from (17) that we can rewrite the JOE regularizer as 

^oeK) = log \a k \ - loge-^l(l - e-(' a -*>W) = e x \a k \ + log \a k \ - log(l - e-^-^OI^I), 

so that for sufficiently large \a k \, log(l — e~^° 2 ~ 0l ^ a ^) « 0, 9\\a k \ » log \a k \, and we have that ^joe ( I a fc I ) ~ 9i\ a k\- 
Thus, for large \a k \, the JOE regularizer behaves like t\ with X f = 2a 2 9\. In terms of the probability model, this 
means that the tails of the JOE mixture behave like a Laplacian with 9 = #i, with the region where this happens 
determined by the value of 9^ — 9\. The fact that the non-convex region of V>joe(-) is confined to a neighborhood 
around could help to avoid falling in bad local minima during the optimization (see Section IV for more details 
on the optimization aspects). Finally, although having Laplacian tails means that the estimated a will be biased 
[16], the sharper peak at allows us to perform a more aggressive thresholding of small values, without excessively 
clipping large coefficients, which leads to the typical over- smoothing of signals recovered using an i\ regularizer. 
See Figure 2 (rightmost column) for an example regularizer based on JOE with parameters 9\ = 20, 62 = 100, and 
the thresholding function it induces. 

The JOE regularizer has two hyper-parameters {9 1^2) which define 6 and that, in principle, need to be tuned. 
One possibility is to choose 9\ and 62 based on the physical properties of the data to be modeled, so that the possible 
values of 9 never fall outside of the range [#1, 62]. For example, in modeling patches from grayscale images with 
a limited dynamic range of [0, 255] in a DCT basis, the maximum variance of the coefficients can never exceed 
128 2 . The same is true for the minimum variance, which is defined by the quantization noise. 

Having said this, in practice it is advantageous to adjust [#i,6y to the data at hand. In this case, although no 
closed form solutions exist for estimating [^1,^2] using MLE or the method of moments, standard optimization 
techniques can be easily applied to obtain them. See Appendix A for details. 
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C. The conditional Jeffreys 

A recent approach to deal with the case when the integral over 6 in the Jeffreys prior is improper, is the conditional 
Jeffreys [ , Chapter 11]. The idea is to construct a proper prior, based on the improper Jeffreys prior and the first 
few no samples of a, (ai, a2, . . . , o n J, and then use it for the remaining data. The key observation is that although 
the normalizing integral J y // I(9)d9 in the Jeffreys prior is improper, the unnormalized prior w{6) = yliW) can 
be used as a measure to weight P(ai, a2, . . . , a no \0), 

w(o) = n«ua*,---^\o)y/m (18) 

/ e p ( a i , a 2 , . . . , a no |0 V ^ (O^C 

It turns out that the integral in (18) usually becomes proper for small no in the order of dim(0). In our case we 
have that for any no > 1, the resulting prior is a Gamma(/^o, A)) distribution with k,o := no and /3o := Y^j=i a j 
(see Appendix A for details). Therefore, using the conditional Jeffreys prior in the mixture leads to a particular 
instance of MOE, which we denote by CMOE (although the functional form is identical to MOE), where the Gamma 
parameters k and f3 are automatically selected from the data. This may explain in part why the Gamma prior 
performs so well in practice, as we will see in Section V. 

Furthermore, we observe that the value of /3 obtained with this approach (/3q) coincides with the one estimated 
using the method of moments for MOE if the k in MOE is fixed to k = + 1 = ^o + 1- Indeed, if computed 
from no samples, the method of moments for MOE gives (3 = (k — l)/ii, with fi\ = ^J2 a j> which gives us 
/3 = n °~^ Q ~ 1 a j = Po-It turns out in practice that the value of n estimated using the method of moments gives 
a value between 2 and 3 for the type of data that we deal with (see Section V), which is just above the minimum 
acceptable value for the CMOE prior to be defined, which is no = 1. This justifies our choice of no = 2 when 
applying CMOE in practice. 

As no becomes large, so does ko = no, and the Gamma prior w(6) obtained with this method converges to a 
Kronecker delta at the mean value of the Gamma distribution, S Ko /p ('). Consequently, when w{9) ~ 8 Ko /f3 (9), the 
mixture f Q P{a\9)w{9)d9 will be close to P(a|fto/A))- Moreover, from the definition of k,o and /3o we have that 
Ko/Po is exactly the MLE of 9 for the Laplacian distribution. Thus, for large no, the conditional Jeffreys method 
approaches the mle Laplacian model. 

Although from a universal coding point of view this is not a problem, for large no the conditional Jeffreys model 
will loose its flexibility to deal with the case when different coefficients in A have different underlying 9. On the 
other hand, a small no can lead to a prior w{6) that is overfitted to the local properties of the first samples, which 
for non- stationary data such as image patches, can be problematic. Ultimately, no defines a trade-off between the 
degree of flexibility and the accuracy of the resulting model. 



IV. Optimization and implementation details 

All of the mixture models discussed so far yield non-convex regularizers, rendering the sparse coding problem non- 
convex in a. It turns out however that these regularizers satisfy certain conditions which make the resulting sparse 
coding optimization well suited to be approximated using a sequence of successive convex sparse coding problems, 
a technique known as Local Linear Approximation (LLA) [ ] (see also [46], [19] for alternative optimization 
techniques for such non-convex sparse coding problems). In a nutshell, suppose we need to obtain an approximate 
solution to 

K 

aj = argmin ||xj — Da^ + A ^(|a/g|), (19) 

k=i 

where ?/;(•) is a non-convex function over R + . At each LLA iteration, we compute by doing a first order 

expansion of ?/>(•) around the K elements of the current estimate ajy , 

= + ^(|ag|) (Vl - |ag|) = ^(|ag|)H + c k , 
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and solving the convex weighted l\ problem that results after discarding the constant terms 

K 

a^ +1) = argmin ||xj - Daf^ + A ^ ^(Kl) 



k=i 

K K 

//H.WlM., I - orrrmlr, _ T> oil 2 _L \ (*) I 



= argmin ||xj - Da^ + A ^^(1^- |)|^fc| = argmin ||xj - Dafj + ^ l a ifc|- ( 20 ) 



a 

/c=l /c=l 



where we have defined A^ := A^'(|ajy |). If ?//(•) is continuous in (0, +oo), and right-continuous and finite at 0, 
then the LLA algorithm converges to a stationary point of (19) [51]. These conditions are met for both the MOE 
and JOE regularizers. Although, for the JOE prior, the derivative ?//(•) is not defined at 0, it converges to the limit 
2(l 2 -0 1 i) W ^ en | ^ | — ^ 0, which is well defined for 62 7^ 0\. If 62 = 0i, the JOE mixing function is a Kronecker 
delta and the prior becomes a Laplacian with parameter 9 = 6\ = 62. Therefore we have that for all of the mixture 
models studied, the LLA method converges to a stationary point. In practice, we have observed that 5 iterations are 
enough to converge. Thus, the cost of sparse coding, with the proposed non-convex regularizers, is at most 5 times 
that of a single l\ sparse coding, and could be less in practice if warm restarts are used to begin each iteration. 
Of course we need a starting point gl^\ and, being a non-convex problem, this choice will influence the approxima- 
tion that we obtain. One reasonable choice, used in this work, is to define ajy = ao, k = 1, . . . , if, j = 1, . . . , N, 
where ao is a scalar so that V>'( a o) = E w [9], that is, so that the first sparse coding corresponds to a Laplacian 
regularizer whose parameter is the average value of 6 as given by the mixing prior w(0). 

Finally, note that although the discussion here has revolved around the Lagrangian formulation to sparse coding 
of (4), this technique is also applicable to the constrained formulation of sparse-coding given by Equation (1) for 
a fixed dictionary D. 

Expected approximation error: Since we are solving a convex approximation to the actual target optimization 
problem, it is of interest to know how good this approximation is in terms of the original cost function. To give an 
idea of this, after an approximate solution a is obtained, we compute the expected value of the difference between 
the true and approximate regularization term values. The expectation is taken, naturally, in terms of the assumed 
distribution of the coefficients in a. Since the regularizers are separable, we can compute the error in a separable 
way as an expectation over each k-th coefficient, C?( a /c) = E Vr ^ q ^(z^) — , where ipk(') is the approximation 
of ipk(') around the final estimate of a&. For the case of q = MOE, the expression obtained is (see Appendix) 



CmoeO/c, K, /3) £^MOE(tt,/3) $k(v) - ^0) logK +(3) + 



a k H 



log/3--. 



a k + /3 

In the MOE case, for n and (3 fixed, the minimum of ( M0E occurs when a k = = ft). We also have 
Cmoe(0) = (^-1)- 1 -/^- 1 . 

The function ( q (-) can be evaluated on each coefficient of A to give an idea of its quality. For example, in 
the experiments from Section V, we obtained an average value of 0.16, which lies between Cmoe(0) = 0.19 and 
min a ^ M0E (tt) = 0.09. Depending on the experiment, this represents 6% to 7% of the total sparse coding cost 
function value, showing the efficiency of the proposed optimization. 

Comments on parameter estimation: All the universal models presented so far, with the exception of the 
conditional Jeffreys, depend on hyper-parameters which in principle should be tuned for optimal performance 
(remember that they do not influence the universality of the model). If tuning is needed, it is important to remember 
that the proposed universal models are intended for reconstruction coefficients of clean data, and thus their hyper- 
parameters should be computed from statistics of clean data, or either by compensating the distortion in the statistics 
caused by noise (see for example [ ]). Finally, note that when D is linearly dependent and rank(D) = R M , the 
coefficients matrix A resulting from an exact reconstruction of X will have many zeroes which are not properly 
explained by any continuous distribution such as a Laplacian. We sidestep this issue by computing the statistics 
only from the non-zero coefficients in A. Dealing properly with the case P(a = 0) > is beyond the scope of 
this work. 
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V. Experimental results 

In the following experiments, the testing data X are 8 x 8 patches drawn from the Pascal VOC2006 testing subset, 7 
which are high quality 640x480 RGB images with 8 bits per channel. For the experiments, we converted the 2600 
images to grayscale by averaging the channels, and scaled the dynamic range to lie in the [0, 1] interval. Similar 
results to those shown here are also obtained for other patch sizes. 

A. Dictionary learning 

For the experiments that follow, unless otherwise stated, we use a "global" overcomplete dictionary D with 
K = 4M = 256 atoms trained on the full VOC2006 training subset using the method described in [35], [36], 
which seeks to minimize the following cost during training, 8 

1 N 

n E { ite - D + A ^-)} + » II dTd IIf > ( 21 ) 

3=1 

II T l|2 

where denotes Frobenius norm. The additional term, /i||D J D|| F , encourages incoherence in the learned 

dictionary, that is, it forces the atoms to be as orthogonal as possible. Dictionaries with lower coherence are well 
known to have several theoretical advantages such as improved ability to recover sparse signals [11], [45], and 
faster and better convergence to the solution of the sparse coding problems (1) and (3) [13]. Furthermore, in [35] it 
was shown that adding incoherence leads to improvements in a variety of sparse modeling applications, including 
the ones discussed below. 

We used MOE as the regularizer in (21), with A = 0.1 and ji = 1, both chosen empirically. See [1], [26], [35] 
for details on the optimization of (3) and (21). 

B. MOE as a prior for sparse coding coefficients 

We begin by comparing the performance of the Laplacian and MOE priors for fitting a single global distribution 
to the whole matrix A. We compute A using (1) with e « and then, following the discussion in Section IV, 
restrict our study to the nonzero elements of A. 

The empirical distribution of A is plotted in Figure 3(a), along with the best fitting Laplacian, MOE, JOE, and 
a particularly good example of the conditional Jeffreys (CMOE) distributions. 9 The MLE for the Laplacian fit is 
9 — N\/ || A|| x = 27.2 (here N\ is the number of nonzero elements in A). For MOE, using (12), we obtained k — 2.8 
and /3 = 0.07. For JOE, 6\ = 2.4 and #2 = 371.4. According to the discussion in Section III-C, we used the value 
k, = 2.8 obtained using the method of moments for MOE as a hint for choosing no = 2 (k,q = no + 1 = 3 « 2.8), 
yielding /3q = 0.07, which coincides with the /3 obtained using the method of moments. As observed in Figure 3(a), 
in all cases the proposed mixture models fit the data better, significantly better for both Gamma-based mixtures, 
MOE and CMOE, and slightly better for JOE. This is further confirmed by the Kullback-Leibler divergence (KLD) 
obtained in each case. Note that JOE fails to significantly improve on the Laplacian mode due to the excessively 
large estimated range #2]- I n this sense, it is clear that the JOE model is very sensitive to its hyper-parameters, 
and a better and more robust estimation would be needed for it to be useful in practice. 

Given these results, hereafter we concentrate on the best case which is the MOE prior (which, as detailed above, 
can be derived from the conditional Jeffreys as well, thus representing both approaches). 

From Figure 1(e) we know that the optimal 9 varies locally across different regions, thus, we expect the mixture 
models to perform well also on a per-atom basis. This is confirmed in Figure 3(b), where we show, for each row 
a fc , k = 1, . . . , K, the difference in KLD between the globally fitted MOE distribution and the best per-atom fitted 
MOE, the globally fitted Laplacian, and the per-atom fitted Laplacians respectively. As can be observed, the KLD 
obtained with the global MOE is significantly smaller than the global Laplacian in all cases, and even the per-atom 
Laplacians in most of the cases. This shows that MOE, with only two parameters (which can be easily estimated, as 

7 http://pascallin.ecs.soton.ac.uk/challenges/VOC/databases.html#VOC2006 

8 While we could have used off-the-shelf dictionaries such as DCT in order to test our universal sparse coding framework, it is important 
to use dictionaries that lead to the state-of-the-art results in order to show the additional potential improvement of our proposed regularizers. 

9 To compute the empirical distribution, we quantized the elements of A uniformly in steps of 2 -8 , which for the amount of data available, 
gives us enough detail and at the same time reliable statistics for all the quantized values. 
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Fig. 3. (a) Empirical distribution of the coefficients in A for image patches (blue dots), best fitting Laplacian (green), MOE (red), CMOE 
(orange) and JOE (yellow) distributions. The Laplacian (KLD=0.17 bits) is clearly not fitting the tails properly, and is not sufficiently peaked 
at zero either. The two models based on a Gamma prior, MOE (KLD= 0.01 bits) and CMOE (KLD= 0.01 bits), provide an almost perfect 
fit. The fitted JOE (KLD=0.14) is the most sharply peaked at 0, but doest not fit the tails as tight as desired. As a reference, the entropy 
of the empirical distribution is H = 3.00 bits, (b) KLD for the best fitting global Laplacian (dark green), per-atom Laplacian (light green), 
global MOE (dark red) and per-atom MOE (light red), relative to the KLD between the globally fitted MOE distribution and the empirical 
distribution. The horizontal axis represents the indexes of each atom, k = 1 , . . . , K , ordered according to the difference in KLD between 
the global MOE and the per-atom Laplacian model. Note how the global MOE outperforms both the global and per-atom Laplacian models 
in all but the first 4 cases, (c) active set recovery accuracy of £± and MOE, as defined in Section V-C, for L = 5 and L = 10, as a function 
of a. The improvement of MOE over i\ is a factor of 5 to 9. (d) PSNR of the recovered sparse signals with respect to the true signals. In 
this case significant improvements can be observed at the high SNR range, specially for highly sparse (L = 5) signals. The performance of 
both methods is practically the same for a > 10. 



detailed in the text), is a much better model than K Laplacians (requiring K critical parameters) fitted specifically to 
the coefficients associated to each atom. Whether these modeling improvements have a practical impact is explored 
in the next experiments. 

C. Recovery of noisy sparse signals 

Here we compare the active set recovery properties of the MOE prior, with those of the £i-based one, on data 
for which the sparsity assumption \Aj\ < L holds exactly for all j. To this end, we obtain sparse approximations 
to each sample Xj using the 4>-based Orthogonal Matching Pursuit algorithm (OMP) on D [ ], and record the 
resulting active sets Aj as ground truth. The data is then contaminated with additive Gaussian noise of variance a 
and the recovery is performed by solving (1) for A with e = CM a 2 and either the t\ or the MOE-based regularizer 
for ip(-). We use C = 1.32, which is a standard value in denoising applications (see for example [27]). 

For each sample j, we measure the error of each method in recovering the active set as the Hamming distance 
between the true and estimated support of the corresponding reconstruction coefficients. The accuracy of the method 
is then given as the percentage of the samples for which this error falls below a certain threshold T. Results are 
shown in Figure 3(c) for L = (5, 10) and T = (2, 4) respectively, for various values of a. Note the very significant 
improvement obtained with the proposed model. 

Given the estimated active set Aj, the estimated clean patch is obtained by projecting Xj onto the subspace 
defined by the atoms that are active according to Aj, using least squares (which is the standard procedure for 
denoising once the active set is determined). We then measure the PSNR of the estimated patches with respect 
to the true ones. The results are shown in Figure 3(d), again for various values of a. As can be observed, the 
MOE-based recovery is significantly better, specially in the high SNR range. Notoriously, the more accurate active 
set recovery of MOE does not seem to improve the denoising performance in this case. However, as we will see 
next, it does make a difference when denoising real life signals, as well as for classification tasks. 

D. Recovery of real signals with simulated noise 

This experiment is an analogue to the previous one, when the data are the original natural image patches (without 
forcing exact sparsity). Since for this case the sparsity assumption is only approximate, and no ground truth is 
available for the active sets, we compare the different methods in terms of their denoising performance. 

A critical strategy in image denoising is the use of overlapping patches, where for each pixel in the image a 
patch is extracted with that pixel as its center. The patches are denoised independently as M-dimensional signals 
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Fig. 4. Sample image denoising results. Top: Barbara, a = 30. Bottom: Boats, a = 40. From left to right: noisy, ^i/OMP, £i/£i, moe/moe. 
The reconstruction obtained with the proposed model is more accurate, as evidenced by a better reconstruction of the texture in Barbara, and 
sharp edges in Boats, and does not produce the artifacts seen in both the £\ and £o reconstructions, which appear as black/white speckles 
all over Barbara, and ringing on the edges in Boats. 
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24.1/26.6 24.2/26.7 
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TABLE I 

Denoising results: in each table, each column shows the denoising performance of a learning+coding combination. 
Results are shown in pairs, where the left number is the psnr between the clean and recovered individual patches, 
and the right number is the psnr between the clean and recovered images. best results are in bold. the proposed 
moe produces better final results over both the £ and £\ ones in all cases, and at patch level for all a > 10. note 
that the average values reported are the psnr of the average mse, and not the psnr average. 



and then recombined into the final denoised images by simple averaging. Although this consistently improves the 
final result in all cases, the improvement is very different depending on the method used to denoise the individual 
patches. Therefore, we now compare the denoising performance of each method at two levels: individual patches 
and final image. 

To denoise each image, the global dictionary described in Section V-A is further adapted to the noisy image 
patches using (21) for a few iterations, and used to encode the noisy patches via (2) with e = CMa 2 . We repeated 
the experiment for two learning variants (l\ and MOE regularizers), and two coding variants ((2) with the regularizer 
used for learning, and £q via OMP. The four variants were applied to the standard images Barbara, Boats, Lena, Man 
and Peppers, and the results summarized in Table I. We show sample results in Figure 4. Although the quantitative 
improvements seen in Table I are small compared to l\, there is a significant improvement at the visual level, as 
can be seen in Figure 4. In all cases the PSNR obtained coincides or surpasses the ones reported in [ ]. 10 

E, Zooming 

As an example of signal recovery in the absence of noise, we took the previous set of images, plus a particularly 
challenging one (Tools), and subsampled them to half each side. We then simulated a zooming effect by upsampling 

10 Note that in [1], the denoised image is finally blended with the noisy image using an empirical weight, providing an extra improvement 
to the final PSNR in some cases. The results in I are already better without this extra step. 
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Fig. 5. Zooming results. Left to right: summary, Tools image, detail of zooming results for the framed region, top to bottom and left to 
right: cubic, £o, £±, MOE. As can be seen, the MOE result is as sharp as £o but produces less artifacts. This is reflected in the O.ldB overall 
improvement obtained with MOE, as seen in the leftmost summary table. 



them and estimating each of the 75% missing pixels (see e.g., [50] and references therein). We use a technique 
similar to the one used in [32]. The image is first interpolated and then deconvoluted using a Wiener filter. The 
deconvoluted image has artifacts that we treat as noise in the reconstruction. However, since there is no real noise, 
we do not perform averaging of the patches, using only the center pixel of Xj to fill in the missing pixel at j. The 
results are summarized in Figure 5, where we again observe that using MOE instead of £q and i\ improves the 
results. 

F. Classification with universal sparse models 

In this section we apply our proposed universal models to a classification problem where each sample Xj is to 
be assigned a class label yj = 1, . . . , c, which serves as an index to the set of possible classes, {Ci, C2, . . . , C c }. We 
follow the procedure of [ ], where the classifier assigns each sample Xj by means of the maximum a posteriori 
criterion (5) with the term — logP(a) corresponding to the assumed prior, and the dictionaries representing each 
class are learned from training samples using (21) with the corresponding regularizer ^(a) = — logP(a). Each 
experiment is repeated for the baseline Laplacian model, implied in the i\ regularizer, and the universal model MOE, 
and the results are then compared. In this case we expect that the more accurate prior model for the coefficients 
will result in an improved likelihood estimation, which in turn should improve the accuracy of the system. 

We begin with a classic texture classification problem, where patches have to be identified as belonging to one 
out of a number of possible textures. In this case we experimented with samples of c = 2 and c = 3 textures drawn 
at random from the Brodatz database, 11 , the ones actually used shown in Figure 6. In each case the experiment 
was repeated 10 times. In each repetition, a dictionary of K = 300 atoms was learned from all 16x16 patches of 
the leftmost halves of each sample texture. We then classified the patches from the rightmost halves of the texture 
samples. For the c = 2 we obtained an average error rate of 5.13% using t\ against 4.12% when using MOE, which 
represents a reduction of 20% in classification error. For c = 3 the average error rate obtained was 13.54% using 
l\ and 11.48% using MOE, which is 15% lower. Thus, using the universal model instead of t\ yields a significant 
improvement in this case (see for example [26] for other results in classification of Brodatz textures). 

The second sample problem presented is the Graz'02 bike detection problem, 12 where each pixel of each testing 
image has to be classified as either background or as part of a bike. In the Graz'02 dataset, each of the pixels can 
belong to one of two classes: bike or background. On each of the training images (which by convention are the 
first 150 even-numbered images), we are given a mask that tells us whether each pixel belongs to a bike or to the 
background. We then train a dictionary for bike patches and another for background patches. Patches that contain 
pixels from both classes are assigned to the class corresponding to the majority of their pixels. 

1 1 http : // w w w. ux . ui s . no/~ tr anden/brodatz . html 
12 http://lear.inrialpes.fr/people/marszalek/data/ig02/ 
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Fig. 6. Textures used in the texture classification example. 




Fig. 7. Classification results. Left to right: precision vs. recall curve, a sample image from the Graz'02 dataset, its ground truth, and the 
corresponding estimated maps obtained with £± and MOE for a fixed threshold. The precision vs. recall curve shows that the mixture model 
gives a better precision in all cases. In the example, the classification obtained with MOE yields less false positives and more true positives 
than the one obtained with £i. 



In Figure 7 we show the precision vs. recall curves obtained with the detection framework when either the l\ 
or the MOE regularizers were used in the system. As can be seen, the MOE-based model outperforms the t\ in this 
classification task as well, giving a better precision for all recall values. 

In the above experiments, the parameters for the l\ prior (A), the MOE model (A M0E ) and the incoherence term 
(ji) were all adjusted by cross validation. The only exception is the MOE parameter (3, which was chosen based on 
the fitting experiment as /3 = 0.07. 

VI. Concluding remarks 

A framework for designing sparse modeling priors was introduced in this work, using tools from universal 
coding, which formalizes sparse coding and modeling from a MDL perspective. The priors obtained lead to models 
with both theoretical and practical advantages over the traditional £o and i\ -based ones. In all derived cases, the 
designed non-convex problems are suitable to be efficiently (approximately) solved via a few iterations of (weighted) 
t\ subproblems. We also showed that these priors are able to fit the empirical distribution of sparse codes of 
image patches significantly better than the traditional IID Laplacian model, and even the non-identically distributed 
independent Laplacian model where a different Laplacian parameter is adjusted to the coefficients associated to each 
atom, thus showing the flexibility and accuracy of these proposed models. The additional flexibility, furthermore, 
comes at a small cost of only 2 parameters that can be easily and efficiently tuned (either (as, 0) in the MOE model, 
or (61,62) in the JOE model), instead of K (dictionary size), as in weighted i\ models. The additional accuracy 
of the proposed models was shown to have significant practical impact in active set recovery of sparse signals, 
image denoising, and classification applications. Compared to the Bayesian approach, we avoid the potential burden 
of solving several sampled sparse problems, or being forced to use a conjugate prior for computational reasons 
(although in our case, a fortiori, the conjugate prior does provide us with a good model). Overall, as demonstrated 
in this paper, the introduction of information theory tools can lead to formally addressing critical aspects of sparse 
modeling. 
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Future work in this direction includes the design of priors that take into account the nonzero mass at a = 
that appears in overcomplete models, and online learning of the model parameters from noisy data, following for 
example the technique in [30]. 
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Appendix 
Derivation of the MOE model 



In this case we have P(a\9) = 9e 0a and w(O\k,0) = ^h)^ 1 /3 K e & e , which, when plugged into (9), gives 



Q(a\/3,K) 
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After the change of variables u :— (a + 0)9 (u(0) — 0, u(oo) — oo), the integral can be written as 
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obtaining Q{a\f3,K) = K,f3 K (a + since the integral on the second line is precisely the definition of 

T(k + 1). The symmetrization is obtained by substituting a by \a\ and dividing the normalization constant by two, 

Q(\a\\P, k) = 0.5^(|a| + 

The mean of the MOE distribution (which is defined only for k > 1) can be easily computed using integration 
by parts, 



(u + fi)^ 1 ) 



du — k/3 



+ 



du 



(3 



In the same way, it is easy to see that the non-central moments of order i are — ^=iy- 

The MLE estimates of n and /3 can be obtained using any nonlinear optimization technique such as Newton 
method, using for example the estimates obtained with the method of moments as a starting point. In practice, 
however, we have not observed any significant improvement in using the MLE estimates over the moments-based 
ones. 



Expected approximation error in cost function 

As mentioned in the optimization section, the LLA approximates the MOE regularizer as a weighted t\. Here 
we develop an expression for the expected error between the true function and the approximate convex one, where 
the expectation is taken (naturally) with respect to the MOE distribution. Given the value of the current iterate 
a® = ao, (assumed positive, since the function and its approximation are symmetric), the approximated regularizer 

is ^W(a) = log(a + p) + ^j+^O - oq). We have 
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13 http://www.di.ens.fr/willow/SPAMS/ 
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Derivation of the constrained Jeffreys (JOE) model 
In the case of the exponential distribution, the Fisher Information Matrix in (15) evaluates to 
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By plugging this result into (14) with 6 = [#i,# 2 ], < 9\ < # 2 < oc we obtain w(9) = ^^y^ \- We now 
derive the (one-sided) JOE probability density function by plugging this w(9) in (9), 



Q(a) = 
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Although Q(a) cannot be evaluated at a = 0, the limit for a — >■ exists and is finite, so we can just define Q(0) 
as this limit, which is 

1 



lim Q(a) 
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Again, if desired, parameter estimation can be done for example using maximum likelihood (via nonlinear 
optimization), or using the method of moments. However, in this case, the method of moments does not provide a 
closed form solution for (9\,9 2 ). The non-central moments of order i are 
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For i = 1, both integrals in (22) are trivially evaluated, yielding m 
can be solved using integration by parts: 
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where the first term in the right hand side of both equations evaluates to for i > 1. Therefore, for i > 1 we obtain 
the recursions nf = ^/i^L l5 ^ = which, combined with the result for i = 1, give the final expression 

for all the moments of order i > 
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In particular, for i = 1 and z = 2 we have 6\ = (ln(# 2 /#i)/ii + # 2 *) 1 > ^2 = (ln(02/0i)/i2 + #i 2 ) 1 ? which, 
when combined, give us 

01 = -i z^ 1 77; x 9 ? 02 = ; ,} /n x 9 • (23) 



-2\-l 



fi 2 - ln(0 2 /0i)Mi 
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One possibility is to solve the nonlinear equation 2 /#i = M2 — in(6> 2 /6>i)^j 
nonlinear equation ^ = ^2 and choosing one of them based on some side information. Another possibility 
is to simply fix the ratio # 2 /#i beforehand and solve for 9\ and 62 using (23). 



Derivation of the conditional Jeffreys (CMOE) model 

The conditional Jeffreys method defines a proper prior w(9) by assuming that no samples from the data to be 
modeled a were already observed. Plugging the Fisher information for the exponential distribution, 1(9) = 9~ 2 , 
into (18) we obtain 

(UT=i Oe- ^- 1 e n -i e -e^ iaj 



w(9) 
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Denoting So = X}j=i a j an< ^ performing the change of variables u := So£ we obtain 

0n o -l e -S o 6 gn O Q no -i e -S o 

W ^ ~ Su no J +o ° u n o-i e - u du ~ rW ' 

where the last equation derives from the definition of the Gamma function, T(no). We see that the resulting prior 
w(9) is a Gamma distribution Gamma(/^o, A)) with = no and /3q = So = X)j=i a j- 
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